Profiling placental DNA methylation associated with maternal SSRI treatment during pregnancy

Selective serotonin reuptake inhibitors (SSRIs) for treatment of prenatal maternal depression have been associated with neonatal neurobehavioral disturbances, though the molecular mechanisms remain poorly understood. In utero exposure to SSRIs may affect DNA methylation (DNAme) in the human placenta, an epigenetic mark that is established during development and is associated with gene expression. Chorionic villus samples from 64 human placentas were profiled with the Illumina MethylationEPIC BeadChip; clinical assessments of maternal mood and SSRI treatment records were collected at multiple time points during pregnancy. Case distribution was 20 SSRI-exposed cases and 44 SSRI non-exposed cases. Maternal depression was defined using a mean maternal Hamilton Depression score > 8 to indicate symptomatic depressed mood (“maternally-depressed”), and we further classified cases into SSRI-exposed, maternally-depressed (n = 14); SSRI-exposed, not maternally-depressed (n = 6); SSRI non-exposed, maternally-depressed (n = 20); and SSRI non-exposed, not maternally-depressed (n = 24). For replication, Illumina 450K DNAme profiles were obtained from 34 additional cases from an independent cohort (n = 17 SSRI-exposed, n = 17 SSRI non-exposed). No CpGs were differentially methylated at FDR < 0.05 comparing SSRI-exposed to non-exposed placentas, in a model adjusted for mean maternal Hamilton Depression score, or in a model restricted to maternally-depressed cases with and without SSRI exposure. However, at a relaxed threshold of FDR < 0.25, five CpGs were differentially methylated (|Δβ| > 0.03) by SSRI exposure status. Four were covered by the replication cohort measured by the 450K array, but none replicated. No CpGs were differentially methylated (FDR < 0.25) comparing maternally depressed to not depressed cases. In sex-stratified analyses for SSRI-exposed versus non-exposed cases (females n = 31; males n = 33), three additional CpGs in females, but none in males, were differentially methylated at the relaxed FDR < 0.25 cut-off. We did not observe large-scale alterations of DNAme in placentas exposed to maternal SSRI treatment, as compared to placentas with no SSRI exposure. We also found no evidence for altered DNAme in maternal depression-exposed versus depression non-exposed placentas. This novel work in a prospectively-recruited cohort with clinician-ascertained SSRI exposure and mood assessments would benefit from future replication.


Discovery cohort
Placental chorionic villi samples were obtained for each of the 64 cases in the discovery cohort. As our previous work has shown that DNA methylation (DNAme) patterns vary between placental sampling sites [1,2], to obtain a representative sample of the whole placenta four independent villi biopsies (1.5-2 cm 3 sites from different cotyledons) were taken from the fetal side of each placenta. DNA was extracted from each of the four villi sites using a standard salting-out procedure [3]. Extracted DNA from each of the four sites per placenta was then combined in equimolar amounts for DNAme analysis.
Combined placental DNA from the 64 placental samples was purified using the Qiagen DNeasy Blood and Tissue kit (Qiagen, CA, USA), and 750 ng of DNA per placenta was bisulfite converted using the Zymo EZ DNA Methylation kit (Zymo Research, USA). Bisulfite-converted DNA was whole-genome amplified, enzymatically fragmented, and hybridized to the array per Illumina's protocol for the Infinium MethylationEPIC BeadChip [4]. The cases from the discovery cohort were assayed as part of a larger 424-sample placental DNAme cohort, with all samples similarly extracted and pooled from four placental sites, as well as purified and bisulfiteconverted before loading on the array chips.
To minimize technical effects of sample processing, the 64 cases from the discovery cohort were run within the same EPIC array batch distributed across 9 EPIC chips. The discovery cohort cases were carefully randomized across the chips and the chip rows (1-8) by SSRI exposure (yes/no) and infant sex assigned at birth (male/female) and were well-distributed by mean maternal Hamilton Depression score (Supplementary Figure 2). The remaining cases from the 424-sample cohort were run in subsequent array batches.
Raw IDAT files from all 424 samples were read directly into R (version 4.1.1) using the minfi and ewastools packages [5][6][7][8]. To exclude the possibility of sample mislabelling, contamination, or misloading, genetic uniqueness of all samples was confirmed by clustering on the 59 explicit genotyping ("rs") probes present on the EPIC array. Placental sex chromosome complement was confirmed to match reported infant sex at birth using the ewastools::check_sex() function [6]. Poor quality probes from the following categories were removed from the dataset: detection p value > 0.01 or beadcount < 3 in > 5% of samples (n=4,351); cross-hybridizing probes and probes with single nucleotide polymorphisms within 5 base pairs of the CpG (n=99,360) [9]; and probes located on the X or Y chromosome (n= 19,637). After sample quality control, probe filtering, and dasen+noob normalization [10], 746,608 autosomal CpGs remained; the 64 discovery cohort samples were then extracted from the 424-sample cohort for separate analysis of SSRI exposure. As a final step prior to analysis, CpGs with non-variable DNAme in the 64 sample discovery cohort samples were removed from the dataset (n=87,572); non-variable DNAme was defined as in [11]: probes with less than a 5% range in DNAme  values between the 10 th -90 th centile that overlapped with CpGs previously reported to be placenta-non-variable in [12]. The final dataset subjected to analysis comprised of 659,036 CpGs in 64 samples.

Replication cohort
Raw IDAT files for 335 samples from GSE75248 were read into R (version 4.1.1) using the minfi and ewastools packages [5][6][7][8]. Genetic uniqueness and sample sex were confirmed to match reported metadata, following the same methods used in the discovery cohort. Poor quality probes were removed from the dataset using the same criteria as applied in the discovery cohort and included: probes with detection p value > 0.01 or beadcount > 3 in > 5% of samples (n=4,312); cross-hybridizing probes and probes with single nucleotide polymorphisms within 5 base pairs of the CpG (n=60,460) [9]; probes located on the X or Y chromosome (n=11,648). After sample quality control and probe filtering, dasen+noob normalization was applied [10], and 411,896 autosomal CpGs remained for analysis. The data were then filtered to only the 34 samples for which gestational SSRI treatment information was available. For replication of differentially methylated CpGs in the discovery cohort, the CpGs to be assessed were selected from the replication cohort dataset, and a linear model was run only on these loci to assess replication of differential DNAme associated with SSRI exposure. This linear model was adjusted for sex, gestational age at birth, PlaNET ancestry, and mode of delivery (vaginal versus Caesarean section); mode of delivery was included as an adjustment covariate as frequency of vaginal delivery differed in association with SSRI exposure in the replication cohort (chi-squared p value 0.016), see main manuscript Table 1. are plotted along the Y axis; points are colored by colored by SSRI exposure (dark yellow = SSRI-exposed, blue = non-SSRI-exposed). Average Δβ across this DMR is +0.06 (higher in SSRI-exposed placentas), and the minimum smoothed FDR calculated by DMRcate across this region was 1.9510 -29 .

Supplementary Figure 5. SSRI exposure is not associated with widespread DNAme patterns in sex-stratified subsets of the discovery cohort. (A) Volcano plot showing the association between DNAme in female placentae and SSRI exposure, and (B) Volcano plot
showing the association between DNAme in male placentae and SSRI exposure. For both plots SSRI refers to selective serotonin reuptake inhibitors, the difference in DNAme (Δβ) is plotted along the X axis and was calculated as Δβ = βSSRI-exposed -βSSRI non-exposed, FDR is shown along the Y axis with more significant (lower FDR) values at the top of the plot. Vertical dashed intercepts demarcate Δβ = 0.03, a horizontal dashed intercept indicates FDR = 0.05. (C) Boxplot of differential DNAme by SSRI exposure status in female samples only (n=31) at cg03905236, in the 5' untranslated region (UTR) of the SH3GL3 gene on chromosome 15. Points are colored by SSRI exposure (blue = SSRI non-exposed, dark yellow = SSRI-exposed) and boxplots indicate mean DNAme β value  one standard deviation. (D) Table of CpGs with differential DNAme at FDR < 0.25 by SSRI exposure status in females only (n=31). Δβ refers to the difference in DNAme between groups and was calculated as Δβ = βSSRI-exposed -βSSRI non-exposed. Pos refers to the genomic coordinates of the CpG in hg19. Figure 6. CpG sites previously reported to have sex-biased DNAme in the human placenta show some evidence for differential DNAme by SSRI exposure status. (A) Scatterplot of DNAme β values  one standard deviation in SSRI-exposed versus SSRI nonexposed samples at cg26136722 in C14orf132 in the discovery cohort; points are colored by SSRI exposure (blue = SSRI non-exposed, dark yellow = SSRI-exposed). (B) Scatterplot of DNAme β values  one standard deviation in SSRI-exposed versus SSRI non-exposed samples at cg22515303 in GTDC1 in the discovery cohort; points are colored by SSRI exposure (blue = SSRI non-exposed, dark yellow = SSRI-exposed). (C) Scatterplot of DNAme β values  one standard deviation in SSRI-exposed versus SSRI non-exposed samples at cg27003571 in GTDC1 in the discovery cohort; points are colored by SSRI exposure (blue = SSRI non-exposed, dark yellow = SSRI-exposed).

Supplementary Figure 7. Testing discovery cohort differentially methylated CpGs for replication in GSE75248. (A)
Boxplot of mean DNAme β values  one standard deviation in SSRI-exposed versus SSRI non-exposed samples at cg14340829 in the replication cohort; points are colored by SSRI exposure (blue = SSRI non-exposed, dark yellow = SSRI-exposed). (B) Boxplot of mean DNAme β values  one standard deviation in SSRI-exposed versus SSRI nonexposed samples at cg20877313 in the replication cohort; points are colored by SSRI exposure (blue = SSRI non-exposed, dark yellow = SSRI-exposed). (C) Plot of DNAme β values in the replication cohort by SSRI exposure status in the SSRI-associated differentially methylated region in the DGKA gene, points are colored by SSRI exposure (blue = SSRI non-exposed, dark yellow = SSRI-exposed) Supplementary Figure 8. Placental epigenetic age acceleration is not associated with SSRI exposure or mean maternal Hamilton Depression score. (A) Placental epigenetic age acceleration (Epi Age Accel) versus mean Hamilton Depression (HamD) score across gestation, colored by SSRI exposure (blue = SSRI non-exposed, dark yellow = SSRI-exposed). (B) Boxplot of placental epigenetic age acceleration separated by SSRI exposure status; boxplots indicate mean  one standard deviation. (C) Scatterplot of mean Hamilton depression score versus placental epigenetic age acceleration with a line of best fit. Figure 9. Mean maternal Hamilton Depression Score is significantly correlated with mean maternal Edinburgh Postnatal Depression Score in the discovery cohort. Mean maternal Hamilton Depression score (HamD) is plotted along the X axis, mean Supplementary Tables   Supplementary Table 1. Results (number of significant CpGs) from SSRI exposure and maternal depression differential DNAme analyses. Model A refers to linear modelling to identify DNAme alterations associated with SSRI exposure in full cohort (n=64), adjusting for mean maternal Hamilton Depression score across gestation. Model B refers to linear modelling to identify DNAme alterations associated with SSRI exposure in depressed samples (n=34, all mean maternal Hamilton Depression score > 8). Model C refers to linear modelling to identify DNAme alterations associated with maternal depression in SSRI non-exposed samples (n=44).